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Abstract 


If quenched to zero temperature, the one-dimensional Ising spin chain undergoes coarsening, 
whereby the density of domain walls decays algebraically in time. We show that this coarsening 
process can be interrupted by exerting a rapidly oscillating periodic field with enough strength to 
compete with the spin-spin interaction. By analyzing correlation functions and the distribution 
of domain lengths both analytically and numerically, we observe nontrivial correlation with more 
than one length scale at the threshold held strength. 
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I. INTRODUCTION 


One of the most important topics in statistical physics is the formation of order. A 
classical nonequilibrium example is provided by the one-dimensional (ID) Glauber-Ising 
model quenched to a zero temperature. It approaches one of the ordered ground states by 


forming larger and larger domains |l|, and this coarsening process has been analyzed in full 
detail (see, e.g., Ref. 2). The relaxation toward equilibrium is very slow: In the absence 
of an external field, the domain walls perform annihilating random walks and the density 

V constrained version 
-6], and such ’’aging” 


accordingly decays as p ~ t 1//2 as time t goes by Q . A kinetical 
also exhibits glassy behavior with anomalously slower coarsening 


can even cease to proceed under steady driving [7]. One may ask if something similar can 
be achieved in the original Glauber-Ising system by driving it with a suitable protocol, 


as has been pointed out in Ref. ]8|. Evidently, a constant external field does not work for 
that purpose because the field will only accelerate the coarsening dynamics by breaking 
the up-down symmetry. If the symmetry is concerned, an alternative protocol would be an 
oscillating field with a short period. This is of particular interest from the perspective of 
interaction of light and matter in the high-frequency regime. The problem becomes highly 
nontrivial, especially when the matter has internal spatiotemporal correlations as a many- 
body system (see, e.g., Refs. 9-ll| and references therein). Due to its ubiquity and often 

the nonequilibrium caused by oscillatory 


|l3, 14]. 


If we 


dramatic consequences as reported in Ref. r], 

driving still remains as an active area of research to be explored further 
consider the ID Ising chain under an oscillating field, one possible scenario is that the 
system has such a large time scale that it simply overlooks the rapid oscillation so that 
the field appears as a small perturbation around the ordered state. Indeed, this has been 
numerically observed in the two-dimensional kinetic Ising model subjected to an oscillating 


field (see, e.g., Ref. ]15|). On the other hand, it also seems plausible that the disordered 
state can remain stable, although energetically unfavorable, just as an inverted pendulum 
is stabilized by fast oscillatory driving Q, [l7]. In this paper, we show that the latter is the 
case when the field strength is greater than or equal to the spin-spin coupling strength. In 
fact, from the dynamic rules defined below, we can readily convince ourselves that all the 
correlations are completely destroyed if the field amplitude exceeds the spin-spin coupling 
strength, which effectively corresponds to an infinite temperature. If the field is weaker 
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than the spin-spin interaction, on the other hand, the up-down symmetry can be broken as 
in an ordered phase because the system cannot escape from the absorbing states with all 
the spins aligned in one direction. Only when the internal and external energy scales are 
equally strong, we observe finite nontrivial correlations and a stationary density of domain 
walls. We will explain this point by calculating correlation functions and the distribution of 
domain lengths both analytically and numerically. 

This paper is organized as follows: An explanation of our model system is given in Sec. [TT] 
Correlation functions and the domain length distribution are analyzed in Sec. IIII[ This is 
followed by a discussion of results and conclusions. 


II. GLAUBER-ISING DYNAMICS 


Let us consider a ID Ising chain with size L under a time-dependent external field H(t). 
The energy function is written as 

L L 




(i) 


i —1 


2—1 


where the spin variable Si can take either +1 or —1 and 4 > 0 is the coupling strength 
between neighboring spins. We will impose a periodic boundary condition by setting S^ + \ = 
S i- The time evolution of this system is assumed to obey the zero-temperature Glauber 


dynamics 


18], which means that every spin flips with the following rate: 


Wi = 


1 if A Ei < 0, 

| if AEi = 0, (2) 

0 if A Ei > 0, 

where A Ei is the energy difference due to a spin flip from Si to —S t . The external field H(t) 
takes a rectangular pulse shape between +H 0 and — H 0 with period 2T, where H 0 > 0 is 
a constant. It is convenient to define r = (t mod 2 T) as a time index within each period. 
Then, the external field is described as follows: 

+H 0 for 0 ^ r < T, 

— H 0 for T ^ r < 2 T. 

As briefly mentioned above, we need to consider competition between the spin-spin interac¬ 
tion and the external driving: If H 0 > J, the field direction solely determines the dynamics, 


m = 


( 3 ) 
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so that the system is equivalent to a collection of non-interacting spins subjected to the 
held. If H 0 < J, on the other hand, the held cannot hip a spin once it is surrounded by 
two other spins in the same direction. As a consequence, the density of domain walls keeps 
decreasing, regardless of the held direction, playing the role of the Lyapunov function in this 
dynamics. This means that the steady states under periodic driving must be the ordered 
ones for H 0 < J, and the deterministic nature of the dynamics suggests that the coarsening 
will not be slower than the held-free case. One can indeed numerically check that the density 
of domain walls decays as p ~ H 1 ^ 2 when T is small but with a smaller prefactor than in 
the absence of H(t). For this reason, we can say that H 0 = J is the most nontrivial point 
due to the interplay between the held and the spin-spin interaction. Henceforth, we will set 
H 0 = J unless otherwise mentioned. 

As is well known, the dynamics can also be analyzed in terms of domain walls. We 
will briehy review three basic processes of the domain-wall dynamics, i.e., pair creation, 
pair annihilation, and propagation, assuming that H(t) = +Hq. First, two domain walls 
are created inside a down-spin domain {• • • 4-4-4, ■ • •} when the held hips the spin in the 
middle with rate |, which results in {• • • 444 Second, the pair-annihilation process 

is possible in two different ways, i.e., {• • • 444 •••} =>■ {• • • 444 •••} with rate 1 or 
{• • ■ 444 ' ’' } =>■ {''' 444 ''' } with rate |. Last, a domain wall propagates when a spin 
hips at a domain boundary, e.g., {• • • 444 • • •} =>• {• • • 444 • • • } with rate 1. In the spin 
language, all these processes tend to align spins along the held direction. Therefore, few 
domain walls exist if the held has been applied for a sufficiently long period. One of our 
primary interests is how the density of domain walls varies in time when the time-dependent 
held in Eq. ([3]) drives the system. 

We can formally describe the Glauber-Ising dynamics by using the transition-matrix 
formulation because it is Markovian. The Ising chain in Eq. ([I]) has N = 2 L microstates. 
Indexing the microstates by a = 1,..., N, we define p a (t) as the probability to find the 
system in state a at time t. The probability distribution can then be denoted as p(t) = 
{pi{t) ,p 2 (t), ...,p N (t)} with a constraint for the conservation of total probability, 54 0: Pa(t) = 
1. One can readily calculate any single-time observable from p(£) in principle, including the 
average domain wall density. The zero-temperature Glauber rates in Eq. ([2]) define an N x N 
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transition matrix M(f) that governs the evolution of p(f) in the following way: 


p(t + At) = M(f)p(£), 


( 4 ) 


where At means a time scale for flipping a single spin. It is reasonable to suppose that every 
spin has a chance to flip during one time step on average, which means that At should be 
proportional to L~ l . The rates are dependent on the external held, so we can distinguish 
the rates under H(t ) = +H 0 from those under — H 0 . It implies that we have to work with 
two transition matrices: 


M (t) = 


M+ if H{t) = +H 0 , 

M" if H{t) = -H 0 , 
which are actually related by a simple coordinate transformation 


( 5 ) 


14]. After one period, 


therefore, the probability distribution at time t = 0 evolves to p(f = 2 T) = M T p(f = 0) 
with My = [(M“) L ] T [(M + ) i ] T . When the system has been entrained by the driving, it 
should be found statistically identical at time t and t + 2 T. This can be regarded as a 
nonequilibrium steady state in a stroboscopic sense. For example, we may observe the 
system at the beginning of every period, i.e., at r = 0 and denote the resulting steady state 
as p oc (r = 0). It is obtained by solving the following equation: 


Poo(t = 0) = M r p 00 (r = 0), (6) 

and the existence of such an eigenvector is guaranteed because both the M + and the M~ are 
stochastic. The steady-state distribution for general r is also obtained in a straightforward 
way. In practice, Eq. (J6|) can be solved only for L < 0(10) because the size of M grows as 
an exponential function of L. From a computational point of view, it is often more efficient 
to sample configurations by using a Monte Carlo method. Figure [T] demonstrates that the 
Monte Carlo sampling precisely reproduces the result from Eq. (J6]) . Our Monte Carlo result 
also shows that the transition-matrix calculation for L = 10 is quite accurate in estimating 
the average density of domain walls in a larger system (Fig. [2]). It implies the following: 
Suppose that we randomly take ten consecutive spins in a large system many times and 
count the frequency of an arbitrary spin configuration i. Our observation suggests that it 
will be more or less similar to pt obtained from the transition-matrix calculation, and it is 
supported by Monte Carlo calculations (not shown). If a large system can be approximated 
as a collection of small ones of L ~ 0(10), it is because the characteristic length scale is 
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FIG. 1. Density of domain walls in the zero-temperature Ising chain of length L = 10, entrained 
by the field in Eq. Q with (a) T = 1, (b) T = 2, (c) T = 5, and (d) T = 10. Each panel shows 
numerically exact results from Eq. ([6|) and Monte Carlo results averaged over 10 5 periods, although 
they are indistinguishable in this plot. 


shorter than 0(10). In other words, this observation suggests weakness of the interaction 
between domain walls. This remark will also be supported by other observations below. 


Another important question in this context is whether a dynamic phase transition (DPT) 
occurs as the half period T is varied. For example, for dimensions higher than one, the 


Glauber-Ising moc 
as T decreases 19 


cl undergoes a symmetry-breaking DPT at a sufficiently low temperature 


20]. Such a DPT is explained by the competition between internal and 


external time scales for relaxation and driving, respectively. However, such a DPT seems 
unlikely in our Ising chain, although the temperature is zero: One clue in Figs. [Hand[2] is that 
the response to +H 0 (0 ^ r < T) is indistinguishable from the one to — H 0 (T ^ r < 2T) for 
any value of T. The magnetization rn = L~ 1 JA S t also oscillates around zero with preserving 
the up-down symmetry for any T (not shown). We will present a more quantitative argument 
for the absence of a DPT by using correlation functions, which we introduce below. 
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FIG. 2. (Color online) Size dependence in the density of domain walls for (a) T = 1 and (b) 
T = 10, obtained by using Monte Carlo calculations, where the solid (dotted) lines represent 
L = 10 (L = 10 3 ). 

III. RESULTS 


A. Correlation functions 

We begin by considering slow driving, e.g., as in Fig. [0(d). One can easily understand 
the behavior of the density of domain walls: At r = 0, for example, the held abruptly 
changes from — H$ to +Hq, whereas most of the spins are pointing downward. The density 
of domain walls thus increases when r is small. As r grows further, however, it is followed 
by a downturn in the density because almost all the spins are aligned in the held direction. 
Then, the held changes to — H 0 again, and all the processes of creation and annihilation of 
domain walls are repeated anew. We will put this description on a more quantitative ground 
by considering correlation functions, and then move on to the case of fast driving. 

Let us recap the time evolution of an individual spin i during At as follows: 

f SAt) with probability 1 — W,At, 

Si(t + At) = < ’ (7) 

I — Si(t ) with probability W,At, 

where W ti is given in Eq. ()2l) . In the limit of At —y 0, the time derivative of magnetization 
and that of the two-point correlation function can be written as 

^§ L = ~ 2 (SiWi), (8) 

and 

d{Si ^ +r) = -2 (SiSi+riWi + Wi+r)) , (9) 
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respectively, where (• • •) means the average over configurations. We now suppose that the 
system experiences +Hq. Enumerating all the possible spin triplets, we can summarize the 
Glauber transition rates in Eq. ((21) as follows: 

W''. = ^[9. + (l-9.)(l-Si)]. (10) 

with gt = -j(l — iSj_i)(l — iSi+i). By substituting Eq. (fTUl) with Eqs. (EJ and (P), we find 
that 

^ = 1(3-2m-ft), (Ha) 

— ,! — — (3 m — 4C r + C r _i + C)+1 — C r - 12 ), (lib) 

at 2 

where m = (Si), C r = (SiS i+r ), and Ci )7 . = (5j_*5*5i +r ) = C r j. Note that we have assumed 
invariance under translation and reflection in the correlation functions. We could also write 
down the evolution of the three-point correlation functions, but it is already obvious that 
the equations will not be closed. To proceed, we need to truncate the endless sequence of 
equations. Our minimalist description is neglecting correlation over a distance greater than 
two, so it reads as 

^ = i(3 -2m-(7,), (12a) 

— .) — — — (3 771 — 4C, + 1 + C 2 — m) , (13b) 

at 2 

^l = ^3m-AC 2 + C 1 ), (12c) 

where we have included the evolution of C 2 , which appears in Eq. (1 11 ai) . This description is 
minimalist in the following sense: Suppose that C- 2 is also neglected. Considering Eq. (II2aj) . 
we see that this makes the evolution of m independent of other correlation functions. Unreg¬ 
ulated by higher-order correlations, it has a fixed point at rn — |, which is unphysical. We 
thus conclude that we need to take into account C 2 at least. Note that our simplified dynam¬ 
ics still admits a fully ordered state with m — C\ — C 2 — 1 as a stationary solution for the 
static field H{t) = +H 0 . Calculating the density of domains walls, p[t) = | [1 — C\(t)\, we 
find a striking agreement between Monte Carlo results and Eq. (TT2T) . numerically integrated 
from an initial condition with m — — 1 and C\ = C 2 = 1 [Fig. [31(a)]. This agreement is also 
consistent with the remark in the previous section that the correlation length is not greater 
than 0(10). 
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FIG. 3. (Color online) Density of domain walls as a function of r. (a) The solid lines are obtained 
by using Monte Carlo simulations for L = 10 3 and T = 15, whereas the dotted lines are obtained 
by numerically integrating Eq. (fT2|). (b) Monte Carlo results for L = 10 3 and T = 1 (solid line) 
remain close to the value ii estimated from Eq. (|12D in the limit of T —> 0 (dotted line). 


Having checked our description for slow driving, we may now consider the opposite limit 
of T —> 0. We assume that this limit restores the up-down symmetry so that m is negligible 
in Eqs. (112bl) and (II 2d) in an average sense. This assumption is supported by the following 
argument: If any remnant magnetization m/0 exists, it means that the system is unable 
to respond to such a rapid held modulation. According to this idea, C 2 would not change 
appreciably upon the held reversal, either. When H(t) = —H 0 , Eq. (I12ap takes a slightly 
different form: 


dm(t ) 
dt 



2 m + C 2 ), 


(13) 


where the right-hand side is written in terms of the same correlation functions based on the 
” freezing” scenario above. The summation of Eqs. (I12al) and (1T3|) expresses the total change 
in m during one period, which must vanish in a steady state. This immediately leads to 
the conclusion that m = 0. In this way, we can argue that a symmetry-breaking DPT, as 
a result of the competition between relaxation and driving time scales, should be absent in 
our system. Put differently, the relaxational time scale does not grow longer than the one 
for driving, and this is consistent with our observation of short correlation lengths. As a side 
remark, we add that the statement of vanishing m contains a subtle point: If m was strictly 
zero all the time, it would imply ^ ^ 0 in Eq. (112 a p or (fl3|) . which is self-contradictory. A 
correct explanation is rather that m will keep changing around zero with a small magnitude. 
In addition to m = 0 in an average sense, the steady-state condition requires that both 
and must also vanish. Solving the set of linear equations resulting from Eqs. (Il2bl) 
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and (j!2cf) . we estimate the stationary density of domain walls as p — \{1 — C\) = b). This 
calculation agrees well with our Monte Carlo result, confirming the existence of domain walls 
in the presence of a fast switching field [Fig. [3(b)]. 


B. Domain statistics 


So far, we have focused on the lowest-order ones among the infinite hierarchy of correlation 
functions, and this turns out to be enough to describe certain average quantities, such as 
the density of domain walls. Now, let us proceed to the detailed statistics of domains to 
gain more information. We begin by considering how domains evolve in time when the field 
is taken to be +Hq. Let P n denote the density of down-spin domains of length n so that 
nP n is equal to the fraction of down spins. We have four mechanisms that affect P n . 


(1) A domain of length n disappears when any of its down spins flips upward. According 
to Eq. ([2]), two spins at the boundary flip with rate 1, whereas the rate is reduced to 
| for the other [n — 2) spins in the bulk. Therefore, the total rate of loss amounts to 
1 x 2 + | x (n — 2) = (| + l), multiplied by P n , for n > 2. Note that this formula does 
not cover the case of a single-spin domain, which disappears via {• • • tit ■ ■ •} =>■ 
{... ttt • • • } with rate 1. 

(2) The density P n increases when a domain of length n+1 shrinks by one at the boundary. 
The contribution is counted as 2 P n+ i because of the two boundary spins. 


(3) We can increase P n by dividing a domain of length l > n + 2 into two pieces in such 
a way that 


s " v ^ 

l 


{••■t|-"ltl“‘IT'"} with rate \ 


2 1 


l—n—1 


{•••tl'"ltl"‘lt'"} with rate 

l—n—l n 


(14) 


If n 7 ^ l — n — 1, this has two different possibilities, each with rate |, so the contribution 
to P n from the domain of length l is equal to Pi. Even if n — l — n — 1, the contribution 
is still Pi because the division creates two domains of length n with rate In total, 
this third mechanism contributes YliLn +2 to P n . 
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(4) The last mechanism is to merge a domain of size l < n — 2 and another with size 
n — l — 1 to create a domain of size n. We can visualize it as 


l n—l—1 


} =>• {• • • || ■ ■ ■ HI • • • 11 • • • } with rate 


(15) 


"N/— 

n 


To evaluate the probability of this event, we need to know the probability of the 
configuration on the left-hand side. The independent-interval approximation (IIA) 
suggests that the lengths can be regarded as totally uncorrelated so that the probability 


can be expressed as 
approximately written as \ Ym=i PiPn-i- 1 - 


21|. The total contribution of this mechanism is thus 


Gathering all these terms, we arrive at 


dP l 

dt 

dPn 

dt 


-Pi + 2P 2 + J2 p h 

1=3 

oo ^ n— 2 

— ^1 + ^ Pn + 2P n +i + P; + - PiP n -i ~i for n > 2. 

l=n +2 1=1 


(16a) 

(16b) 


The next step is to consider the dynamics of up-spin domains with keeping the same field 
direction. Similar to P n , we define Q n as the density of up-spin domains of length n. We 
have five mechanisms to affect Q n . 


(i) In the first mechanism, a domain of a single up spin evaporates via {• • • ||| • • • } 

{• • • HI • • •} with rate |. This takes place only for n = 1. 

(ii) Again, this second mechanism applies only to n — 1. A domain of length 1 can be 

created via {• • • ||| •••} {• • • ||| •••} with rate |. For this to happen, we 

have to pick up a down spin surrounded by two other down spins. For a down-spin 
domain of size l, we have l — 2 such spins. Therefore, we compute this contribution as 
| 5 (^ 3 (Z — 2 )Pi. Note that the dynamics of Q n is coupled to that of P n . 

(iii) The third mechanism describes a loss due to the growth from length n. The domain 
can grow to the left or right, each with rate 1 , so the contribution becomes — 2 Q n . 

(iv) A domain of length n can be gained from the growth process as well, when a domain 
of length n — 1 expands to n by flipping a spin upward at the boundary with rate 1 . 
However, we cannot simply write it as 2Q n _i because the spin flip may merge this 
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domain with another. For example, if we look at the left boundary, the following 
process creates a domain of length n: 


71—1 

whereas the following does not: 


(17) 


• 1U t ■ ■ • t-1- ■ ■ •} =1 {• ■ ■ t tt ■ 14- ■ ■ 


(18) 


n— 1 


In short, it depends on the direction of the spin drawn as a double arrow on the 
leftmost side. In a similar spirit to the IIA, we assume a well-defined probability d> 
for the spin to point downward so that the contribution becomes 2Q n _ 1 <3>. 


(v) The last mechanism is to merge two up-spin domains, one with size l and the other 
with size n — l — 1 as follows: 


{• • • 11 ''' 111 '' • t + ''' > =* {'' ~ 1 1 ''' ttt''' 1 4- ''' } with rate 1. (19) 

l n—l—l n 

As before, we resort to the IIA to estimate the contribution as Y^i=i QiQn-i-i- 


To sum up, we have derived equations for Q n as 


dQi 

dt 


1 1 

_ -jQi - 2 Qi +2 - 2 )Pu 


dQn 

dt 


n —2 

—2Q n + 2Q ri _i$ + QiQ n -i -1 

i=i 


for n > 2. 


(20a) 

(20b) 


Even if H{t) = —H 0 , we can derive essentially the same as Eqs. (TT6|) and (l20]h provided 
that the variable Q n indicates domains in the direction of the field, whereas P n does in the 
opposite direction. 

Suppose that T is so short that the down-spin domains are effectively subjected to both 
Eqs. (TT6T) and (1201) . The steady-state condition implies that ^ = 0 for every n > 1. 

As the up-down symmetry is restored, we may also equate every Q n with P n with setting 
$ = |. We finally end up with the following set of equations: 

0 =--Pi + 2P 2 + ^ y-P/, (21a) 

1=3 z 

oo g n— 2 

0 = P n -1 — (3 + —^ P n + 2P n+ i + Pi + — PiP n -i -1 for n > 2. (21b) 

l=n+2 1=1 
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FIG. 4. (Color online) Domain size distribution, (a) The dotted line shows the trial solution 
P n = Az n with A ~ 0.236 629 and z — 0.615 633 that approximately solves Eq. (|2T|) . The crosses 
represent a numerically exact solution of Eq. (1211) truncated at n = 50. The circles are Monte Carlo 
results for the Ising chain of length L = 10 . The inset shows the same data in a semi-logarithmic 
plot, (b) Correction from the simple exponential form [Eq. (1221) ]. We estimate B ~ 0.01 and 
w — 0.75 by fitting the data on a logarithmic scale. 


Our trial solution is an exponential distribution, i.e., P n = Az n with positive constants A 
and z < 1. Substituting this into Eq. (121 aft , we obtain z — 0.615 633. It is worth noting that 
z would be equal to | if all the correlations were destroyed as in the infinite-temperature 
limit. Although this trial solution does not exactly solve Eq. (I21bj) . we can estimate the 
amplitude A ~ 0.236 629 by taking n —> oo. As a cross-check, we truncate Eq. (12T1) by 
setting P n = 0 for n > 50, and solve the 50 coupled equations for Pi,..., P 50 simultaneously. 
It confirms the validity of our trial solution even for small values of n as shown in Fig. 0](a). 
Of course, we have to ask ourselves whether Eq. (12 If) . involved with several uncontrolled 
approximations, correctly describes the domain dynamics. This is checked by simulating 
an Ising chain of length L = 10 4 to sample the domain length distribution. As depicted in 
Fig. H](a), the result shows that Eq. (l2Tj) works qualitatively but tends to underestimate P n 
when n is large. The correction from Az 11 reveals another length scale in the following form: 


P n - Az n ~ Bw 11 , (22) 

with B ~ 0.01 and w ~ 0.75 [Fig. 01(b)] . The second length scale corresponds to roughly 
four lattice spacings, about twice larger than the first one, but its the origin is not fully 
understood yet. 
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IV. SUMMARY 


To summarize, we have considered the zero-temperature Glauber dynamics in the ID Ising 
chain driven by rectangular pulses of period 2 T and strength equal to J. We have argued 
that the driving interrupts the coarsening so that the density of domain walls converges to 
a nonzero stationary value ~ | in the limit of fast driving. We have also calculated the 
steady-state distribution of domain lengths in the same limit by using the IIA, and the 
result indicates the existence of finite nontrivial correlation. Moreover, our Monte Carlo 
calculation shows that the actual density is higher than expected from simple exponential 
decay, revealing the existence of the second length scale, about twice larger than the first 
one. 
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